A forest-fire analogy to explain the fo-value of the Gutenberg-Richter law for 

earthquakes 
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The Drossel-Schwabl model of forest fires can be interpreted in a coarse grained sense as a model for 
the stress distribution in a single planar fault. Fires in the model are then translated to earthquakes. 
I show that when a second class of trees that propagate fire only after some finite time is introduced 
in the model, secondary fires (analogous to aftershocks) are generated, and the statistics of events 
becomes quantitatively compatible with the Gutenberg Richter law for earthquakes, with a realistic 
value of the b exponent. The change in exponent is analytically demonstrated in a simplified 
percolation scenario. Experimental consequences of the proposed mechanism are indicated. 

PACS numbers: 



Many physical systems react upon a continuous input 
of energy, by releasing the accumulated energy in dis- 
continuous bursts, that are called avalanches in general, 
occurring when some threshold condition is reached. Ex- 
amples range from sand piles, to magnetic domain inver- 
sions in ferromagnets, stress release on the earth crust in 
the form of earthquakes, and many others. Many times 
these avalanches display a broad size distribution, typi- 
cally following a power law, that is taken as a manifes- 
tation of the lack of intrinsic spatial scale in the system. 
Although in real systems this fact can only be approx- 
imate, it is usually taken as the ideal target for theo- 
retical models that try to mimic the observed behavior. 
This concept has gone under the name of self-organized 
criticality 1 . However, in order to describe a natural phe- 
nomenon that displays power law distributions, it is not 
necessary for theoretical models to be self organized crit- 
ical. To be considered realistic, they need only to display 
a phenomenology consistent with observation, which is 
always limited, in spatial and temporal scales. 

For the case of earthquakes (EQs), the experimen- 
tal evidence indicates that they follow the so called 
Gutenberg-Richter (GR) law^, stating that the num- 
ber of EQs in a given magnitude interval dM follows 
the empirical relation N(M) ~ 10~ bAI . This has been 
observed to be accurate in a large range of magnitudes 
with values of b ranging from 0.8 to 1.2. When expressed 
in terms of the seismic moment S (such that, up to a 
constant, M = |log 10 S'), this law can be recast in the 
form N(S) ~ S~ T with r = 1 + §&, so realistic values of 
t are in the range 1.5-1.8. Statistical models devised to 
describe EQs 4 (most remarkably, the Burridge-Knopoff 5 , 
and Olami-Feder-ChristensenS. models) typically consider 
the case of a single planar fault. GR law has been found 
in this kind models, and even a correct value of r has been 
obtained in particular cases, typically by adjusting some 
parameters of the model, what goes against the robust- 
ness observed experimentally in very different conditions. 

In the last years it was realized that even if a correct 
GR law is reproduced, the sequences of events generated 
by the models mentioned above are unnatural in the fol- 



lowing sense. Real EQ sequences display the peculiar 
feature of aftershocks (ASs), namely events that appear 
exclusively as a consequence of previous large events, and 
that are not directly correlated to the driving in the sys- 
tem. The existence of ASs indicates that there is some 
other important time scale at play, in addition to the one 
imposed by the external driving. The physical mecha- 
nisms associated to this internal time scale is related to 
plastic effects, creep dynamics, etc., and I refer to them 
in general as relaxation mechanisms. 

Models including this kind of relaxation have recently 
been introduced^— . It was shown that they are able to 
produce ASs with realistic features. Most remarkably, 
when ASs occur, the exponent of the GR law was ob- 
served to change, and set to a value close to r ~ 1.7, 
without fitting any parameter. The reason of this fact 
has not been adequately discussed up to now. 

An appropriate scenario for this discussion turns out 
to be the Drossel-Schwabl model of forest fires 1 ^. Al- 
though introduced to model a very different problem, it 
is recognized that a coarse grained version of the model 
may qualitatively represent the evolution of the stress in 
a single planar fault system.— ~— I summarize here the 
rules of the DS model for completeness. In an initially 
empty, two dimensional square lattice, sites are occupied 
randomly by trees, one every time unit to- After r tree 
insertions, a lightning event occurs at some random posi- 
tion, burning the tree that is located there, and all trees 
that can be reached from this site through nearest neigh- 
bors occupied sites. This instantaneous burning defines 
the fires in the system, their size (the number of burnt 
trees) being denoted by S. In the limit in which r — > oo 
there are fires with arbitrarily large values of S. There 
is evidence that the system does not become truly criti- 
cal in this limi t 14 ! 15 , however a size distribution of fires 
with an exponent around r ~ 1.2 has been consistently 
observed for a wide range of sizes, and this is sufficient 
for the analysis presented here. 

A coarse grained variable in the model is the density 
of trees over some finite size region. This density can 
take values from to 1, and can be thought to represent 
the stress state of the plate. The random addition of 



trees can be considered, in the coarse grained sense, as 
a smooth increase of the stress in the system, caused by 
external tectonic loading. When stress is high in some 
portion of the system, it can be abruptly reduced by a 
random lightning to some small value, through the occur- 
rence of a fire that is interpreted in this view as an EQ. 
The microscopic stochastic laws of the DS model prevent 
it for reaching a globally synchronized oscillatory state. 
The wide distribution of fires in the DS model corre- 
sponds to a wide distribution of earthquakes. However 
the exponent r ~ 1.2 observed in the DS model does not 
correspond to the exponent for EQs, around 1.7. Also, 
as the system lacks an internal dynamics, ASs are not 
produced. However, a physically motivated mechanism 
to produce ASs also produces a realistic GR law, as I will 
explain now. 

I propose a modified DS model by introducing a sec- 
ond species of trees, called B trees (the original ones are 
noted A trees). The rules of the modified model are the 
following. In the insertion step, each time a site is chosen 
for insertion and if this site is empty, it becomes occu- 
pied by an B tree with probability w, and by a A tree 
with probability 1 — w. The value of w will be assumed 
to be small (w <C 1) . In the burning step, A trees burn 
and propagate fire instantaneously to all neighbors, irre- 
spective if they are A or B trees. However, B trees are 
supposed to propagate fire to neighbors only after some 
random time of the order of a new time scale t\. This 
time scale t\ is supposed to be much smaller than the 
insertion time ioj yet much larger than the propagation 
of fire through A trees. An example of the time sequence 
of fires in this modified DS model is shown in Fig. [TJ In 
the initial configuration, there is a given distribution of 
A and B trees partially filling the lattice, and a lightning 
strikes at some random site. This produces an instanta- 
neous fire that burns all trees connected to the original 
site, and that "activates" the B trees highlighted in (b), 
which retard the propagation for some time of the order 
of t\. As these times for individual B trees are considered 
to be random, for the implementation I simply choose one 
of the active B trees randomly and start a secondary fire 
(c) . The process is continued until no tree is burning (d) . 
At this stage, the insertion process is continued. 

Time sequences of fires consist now of "clusters" of 
events that are triggered by lightnings, and are separated 
by the time to . Each cluster is formed by events separated 
by time intervals of the order of t\. The precise values 
of these times will be important in a realistic considera- 
tion of ASs. Here I simply consider that these secondary 
fires are exactly triggered one every time t\. In the EQ 
analogy, secondary fires are the ASs, and active B trees 
are their epicenters. The typical number of ASs observed 
depends on the parameter w, fixing the ratio between B 
and A trees at insertion. Clusters of events correspond to 
the initial shock and all its ASs. The assumed condition 
ti/to — > allows to clearly identify ASs in the model, a 
situation that is not fulfilled in actual seismicity. 

A temporal sequence of events for the modified DS 
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FIG. 1: (Color online) Propagation of fires in the modified DS 
model. White (gray) circles indicate A (B) trees, (a) In some 
initial configuration, a lightning strikes the site indicated by 
the arrow, (b) The state after the instantaneous burning of 
the connected cluster of A sites. The two B trees that are 
highlighted as black have become "active", (c) One of the 
active B sites is chosen at random and propagates fire to its 
neighbor trees. Additional B sites may become active, (d) 
Final stage when there are no more active B sites. Along 
the whole process, fires of sizes 5, 4, 4, and 6 have occurred 
(contours are highlighted in (d)). The total size of the burnt 
cluster is thus 19. 

model is presented in Fig. [2] The clustered structure 
of the sequence is apparent. Secondary fires to a given 
initial lightning occur at the same value oit/to, and thus 
they appear on the same vertical line. The plot in panel 
(b) shows the same events, but plotted as a function of 
the internal time within the cluster t/t±, with the value 
of t set to at the initial lightning of every cluster. Size 
distribution of events are presented in Fig. |3]with open 
symbols. The presence of B trees generates a size dis- 
tribution with a t exponent in the range 1.7-1.8, well 
different from the original slope (r ~ 1.2) (more system- 
atic results are presented in the supplementary material) , 
and comparable to the value of the actual GR law. The 
slope change generated by the presence of B trees (by the 
existence of ASs, in the seismic perspective) is the main 
results of this work, and I will now explain its origin. 

First of all, I notice that in addition to calculate the 
size distribution of individual events, it is also interesting 
to calculate the size distribution of clusters. Namely, I 
define the cluster size Sc as the sum of sizes of the initial 
event and all its secondary fires. The size of the cluster 
Sc is exactly what we had obtained as a single event if B 
trees also propagate fire instantaneously, or equivalently, 
if there are no B trees at all (w — 0). In fact, the dis- 
tribution of clusters, plotted in Fig. [3] with full symbols, 
coincides with the distribution in the original model 14 . 
Thus we see that the effect of B trees is to fragment 
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1000 2000 t/t 3000 

FIG. 2: (a) Temporal sequence of fires, in the modified DS 
model on a system of size 4000 x 4000 system, with r = 5x I0 4 
and w = 10 -3 . Events that occur at the same value of t/to 
form what I have called a cluster, corresponding to the initial 
fire, generated by lightning, and all secondary fires propagated 
through B trees. In (b), the same events are plotted as a 
function of its internal time within the cluster. This time is 
set to for the first event of the cluster, and increases by ti, 
at each secondary fire. 



the clusters in pieces that burn at different times. This 
fragmentation process, at the same time that generates 
secondary fires, produces the change in the exponent of 
the size distribution. 

It is useful to have an unambiguous characterization 
of this effect in a simpler model than the DS forest fire. 
In fact, in the DS model the spatial distribution of trees 
is highly correlated, and rigorous results for the size dis- 
tribution of avalanches are not available. However the 
analysis becomes much simpler if we consider the propa- 
gation of fires on a totally uncorrelated distribution of A 
and B trees, with spatial densities pa and ps (which now 
have to be fixed by hand). This problem can be studied 
with the tools of standard percolation theory*^ If no B 
trees are present (ps = 0), events correspond to those 
that occur in a site percolation problem with probability 
Pa- As pa — > p c — 0.5927, the events become power law 
distributed with an exponent 187/91 - 1 ~ 1.055.— In 
the presence of B trees, we can calculate again the dis- 
tribution of single events and the distribution of clusters. 
Clusters are distributed as if B trees are not present, 
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FIG. 3: (color online) Statistics of events for simulations with 
r = 10 3 , 10 4 , 10 5 , and 10 6 (from left to right), and w = 0.01 
(system size is up to 20000 x 20000). Lower curves correspond 
to individual events, whereas the upper ones are the results 
for the statistics of clusters. Statistics of clusters coincides 
with the result for single events in the case w = 1. 

i.e, they correspond to normal percolation, with a cut 
off fixed by the total density pa + Pb (< Pc, to avoid 
an infinite fire). For single events instead, we have the 
same fragmentation effect discussed in the context of the 
DS model. The clusters become diluted with a fraction 
Pb/(pa + Pb) of B trees, generating a fragmentation ef- 
fect. The size distribution of single events is the size dis- 
tribution of these fragments. The effect of fragmenting 
a percolation cluster by the removal of a small fraction 
of sites was studied ir>i£, where it was shown that close 
to the percolation threshold, the removal of a single site 
generates fragments that are distributed according to a 
power law with an exponent <p, the fragmentation expo- 
nent, that it was found to be (f> ~ 1.60 in two dimensions. 

Results of numerical simulations in the two species ver- 
sion of the percolation problem (see results in the sup- 
plementary material) confirm this view: the inclusion of 
a small fraction of B trees changes the distribution of 
events to a new power law with a r ~ 1.60. A per- 
fect power law is obtained in the limit pa + Pb Pc, 
Pb/ [Pc — (pa + Pb)} — > const. If ps is kept finite as 
Pa + Pb Pc, a non critical distribution is obtained, 
with maximum size of the fragments controlled by the 
value of pa ■ 

The present analysis of the modified DS model and the 
percolation limit points out a simple mechanism that can 
justify the observed value of the b exponent of the GR 
law in the context of a single planar fault situation. It 
suggests that the observation of this value is closely re- 
lated to the presence of aftershocks. It also points out a 
few possibilities to observe experimental signatures that 
would be a consequence of this mechanism. One of them 
is the fact that the statistics of clusters must display a 
less steep decay than that of individual EQs. An ac- 
tual verification of this fact would imply the unambigu- 
ous classification of events in clusters, which is a non 
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trivial task, since the to and t\ time scales are not well 
separated in actual seismicity. However, there have been 
promising attempts to group the EQs in clusters accord- 
ing to appropriate metrics in the space-time-magnitudc 
parameter space (see for instance^), that may serve to 
this purpose. A careful study of experimental data is 
needed to advance in this direction. 
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FIG. 4: (color online) (a) Partial size distribution of events 
in the modified DS model (r = 10 5 , w — 0.01), restricted 
to the indicated events in the clusters. The distribution has 
an exponent close to that of the original model for the first 
event, and becomes progressively steeper as successive sec- 
ondary fires are considered, (b) An actual example of this 
phenomenon (details are provided in the supplementary in- 
formation) . The lower points display the distribution of EQs 
in Southern California during a 20 years time period. The 
upper curves correspond to the distributions of events that 
occurred in the two time windows indicated, after EQs of 
magnitude 4.5 or larger. The behavior is qualitatively similar 
to that in (a). Straight lines are for reference only. 

A second signature indicating that the proposed mech- 



anism may be relevant in actual seismicity is the follow- 
ing. If we take the events as a function of its internal 
time within an AS sequence (Fig. [U(b)) and calculate 
the size distributions in small windows of the internal 
time, we obtain the results in Fig. |4j It is seen that the 
decaying of the distribution has an exponent similar to 
that of the whole clusters for the first event, and becomes 
progressively steeper as successive events within the se- 
quence are considered. This behavior can be understood 
by considering again the simpler percolation scenario. In 
that case, the very first event in each cluster is taken 
from a distribution with a density pa of occupied sites, 
and this has to follow a decaying law with the smaller r 
(~ 1.055 in the percolation case). Successive ASs pro- 
gressively drive the distribution towards the new, larger 
r value. 

This effect is suitable for experimental verification. In 
fact, EQ catalogs are dominated, after large magnitude 
events, by the ASs to those events. Although other events 
that are not ASs certainly occur, they are few in the first 
hours after large events. So they make a minor contri- 
bution that can be neglected in a first approach. In Fig. 
E^b) I show a preliminary analysis of the seismicity in 
Southern California over a 20 years period. The ten- 
dency to display a weaker decay right after large events 
is clear from this figure, and gives support to the present 
proposed mechanism. 

Summarizing, I have considered an analogy between a 
forest fire model and the stress in a single planar fault. 
The inclusion of a second tree species that delays prop- 
agation of fire, was shown to be analogous to include 
the possibility of ASs in the seismic counterpart. In the 
modified model, EQs (or fires) appear in clusters formed 
by an initial event and all its ASs. The change of expo- 
nent of the GR law is a consequence of the fact that the 
clusters are fragmented by the existence of ASs. An ap- 
proximate case that can be analyzed with standard per- 
colation theory gives analytical support to this scenario. 
Also, observable signatures of this mechanism have been 
pointed out, and a preliminary analysis of actual seismic 
data where one of these signatures shows up has been 
presented. 
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Supplemental material to: A forest-fire analogy to explain the b-value of the Gutenberg-Richter law 
for earthquakes 

A: Size distribution of events for an uncorrelated distribution of A and B trees 

Let us assume a random distribution of A and B trees on a square lattice, with densities pa and ps- An equivalent, 
more convenient set of variables is: p = pa + Pb, and e = pb/pa (note that s is related to w introduced in the main 
paper by w = e/(l + e)). We want to calculate N(S,p,e), the probability distribution of fires as a function of their 
size S 1 , for given values of p and e. This function displays scaling properties as p — > p c ~ 0.5927 and for small e. An 
appropriate scaling form is given by 

N(S,p,e) = (S(p c -p)V°, £ /( Pc , (1) 

where a parameter dependent normalization factor has not been included. The value of the exponent a is exactly 
known in 2D percolation, and is given by a — 36/91. Numerical results adjust very well to this form. They are shown 
in Fig. 

The numerical simulation in the present percolation case, is made in the same way as described in Fig. 1 for the 

modified DS model, but using an uncorrelated distribution of A and B trees. In order to to speed up the simulations, 

only the trees that are burnt are actually generated. This is done by starting a cluster at the origin, and sequentially 

extending it, introducing A and B trees with their corresponding probabilities pa and ps- 

The distribution obtained for ps = corresponds to the usual site percolation case. In particular, N°(x, 0) ~ 
2,-96/91 f or sma n x<j so f or p _ i p c we nave 

N(S,p,0) = N o (S(p c -pf,0) ~ S- 96 ' 91 (2) 

As ps becomes different from zero, the size distribution of events departs from the result of standard percolation and 
starts to develop a region at low events sizes where a distribution according to the fragmentation exponent <f> ~ 1.60 
is displayed, i.e, for constant e, and for p — > p c we get 

N(S,p,e)~S-* (3) 

and this behavior remains correct for S < (p c — p)~° ■ This means that a strict power law with the fragmentation 
exponent is obtained in the double limit p — > p c , s/(p c — p) —> const. ^ 0, i.e., the density of B trees has to be 
reduced to zero at the same time that the total density goes to p c . If the limit p — > p c is taken keeping a finite value 
of pb the limiting distribution has an exponential cutoff at large events sizes that is controlled by the value of pa < Pc- 

B: Simulations of the modified DS model 
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FIG. 5: Results of simulations in the percolation case, at two different values of total density p = pa + pa , and for different 
relative density of A and B trees e = pa/ pA- The results follow the scaling form presented in Eq. ([TJ (with a — 36/91), which 
allows to conclude that a pure power law with an exponent <f> ~ 1.60 is obtained in the double limit p — > p c , e/(p c — p) — ► 
const. =fc 

Here I present additional results for the modified DS model. They are shown in Fig. S£j and correspond to 
simulations with progressively larger values of r, and different values of the probability of B trees at insertiuon w. For 
non-zero w, the main dependence of N with S is absorbed by plotting NS T and using r = 1.8. A tentative scaling is 
made by plotting the horizontal axis as S/r x , using A = 1.08 which is the expected scaling in the original DS model. 
We see a non-perfect scaling of the results, with systematic deviations that remain even for the largest values of r 
used, as it was the case of the original DS mode l 14 ' 15 . 

For w in the approximate range 10 -3 — 10~ 2 there is a progressively wider region when r is increased in which 
NS 1 ' 8, is practically constant, meaning that a distribution with the modified exponent N(S) ~ S* -18 is obtained. 
By comparing Figs. 33 and S£]we can conclude that, in addition to the lack of true criticallity, a further difference 
between the results for the modified DS model, and those in the percolation limit is that the more accurate scaling 
of the data in Fig. £|6]is obtained by keeping a fixed value of w as r is increased, contrary to the percolation case, in 
which the density of B trees has to be reduced to zero at the same time that p —> p c . The origin of this difference 
(and whether it persists for larger values of r or not) deserves further investigation, but it can be anticipated that 
the answer is contained in the structure of the fires in the DS model compared with the normal percolation clusters. 
More systematic simulations in larger systems would be helpful to adress this issue, but I note that the distribution 
of events following N(S) ~ £ -1 " 8 extends already through a range of S values that is comparable to the experimental 
range for EQs. 

C: Details on the data in Fig. 4(b) 

Here I give details on the data analysis to produce Fig. 4(b). The starting point is the list of all EQs (chara- 
terized by magnitude and time) in the region of Southern California (between 115° - 119° W, and 32° - 37° N), 
that occurred between 1/1/1980 an d 5/20/2008 with magnitude larger that 2.5 (according to the ANSS database, 
|http:/ /earthqu ake.usgs !gov/monitoring/anss/| about 26000 events are counted). The number of EQs in small magni- 
tude windows were calculated and the results are shown by the squares. Next, the 108 events with magnitude larger 
that 4.5 were considered, and the statistics was limited to the EQs with time of ocurrence between tA and ts after 
any of these 108 large EQs. The results are plotted as circles (tA = 0, ts — 30 min, 970 events are counted) and 
triangles (tA = 30 min, tg = 90 min, 930 events are counted). Note that these two partial distribution are limited to 
magnitudes lower than 4.5, since larger events are counted as main shocks. 
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FIG. 6: Results of simulations of the modified DS model as a function of r, and for different values of w. Groups of curves at 
a same value of w have been displaced vertically by the same factor, for clarity. 



